home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 January / macformat-020.iso / Shareware City / Education / RLaB / toolbox / ode4.r < prev    next >
Encoding:
Text File  |  1994-07-18  |  2.6 KB  |  141 lines  |  [TEXT/ttxt]
  1. //
  2. //  This is an adaptive 4th-order Runge-Kutta integrator.
  3. //  See Numerical Recipes.
  4. //
  5.  
  6. static (rk4, collapse);
  7.  
  8. ode4 = function( ydot, t0, tf, y0, output, tol, h, yscale ) 
  9. {
  10.   local(t, y, dydt, scale, hh, yt, ys, yout, I, ...
  11.           tout, cnt, output_max, TINY, errmax );
  12.  
  13.   TINY = 1e-4;
  14.  
  15.   //  Supply defaults.
  16.   if (!exist (tol)) { tol = 1e-3; }
  17.   if (!exist (h)) { h = abs( tf - t0 ) / 1000.0; }
  18.  
  19.   //  Initialize.
  20.   t = t0;
  21.   h = h * ( tf - t0 ) ./ abs( tf - t0 );
  22.   y = y0[:];
  23.   I = 2;
  24.  
  25.   // Save the initial  conditions
  26.  
  27.   if (!exist (output)) 
  28.   {
  29.     yout.[1] = [t, y.'];
  30.   else
  31.     yout.[1] = output(t,y).';
  32.   }
  33.  
  34.   while ( (t-tf)*(tf-t0) < 0.0 ) 
  35.   {
  36.  
  37.     dydt = ydot( t, y );
  38.     if (!exist (yscale)) { scale = abs(y) + abs(h*dydt) + TINY; }
  39.  
  40.     //    If the step can overshoot end, cut the stepsize.
  41.  
  42.     if ( (t+h-tf)*(t+h-t0) > 0.0 ) { h = tf - t; }
  43.  
  44.     //    Take a step
  45.     ys = y;
  46.  
  47.     while ( 1 ) 
  48.     {
  49.       //  Take two half steps.
  50.  
  51.       hh = 0.5 * h;
  52.       yt = rk4( ydot, t, hh, ys, dydt );
  53.       y =  rk4( ydot, t+hh, hh, yt, ydot( t+hh, yt ) );
  54.  
  55.       //  Check for step too small.
  56.  
  57.       if ( t+h == t ) 
  58.       {
  59.           fprintf("stderr", "rlab: run time error: Step size too small in ode4\n");
  60.           yt = 0.0;
  61.         errmax = 1.0;
  62.         tf = t;
  63.       else
  64.     //  Take the full step.
  65.           yt = y - rk4( ydot, t, h, ys, dydt );
  66.     
  67.     //  Evaluate accuracy.
  68.         errmax = norm( abs( yt ) ./ scale, "I" ) / tol;
  69.       }
  70.  
  71.       //  Success?
  72.       if ( errmax <= 1.0 ) { break; }
  73.  
  74.       //  Try again.  Reduce step size, but by no more than
  75.       //  a factor of about 50.
  76.  
  77.       if ( errmax > 4e6 ) 
  78.       {
  79.           h = 0.02*h;
  80.       else
  81.         h = 0.9*h*errmax.^-0.25;
  82.       }
  83.     }
  84.  
  85.     t = t + h;
  86.  
  87.     //    Take care of 5th order.
  88.     y = y + yt./15.0;
  89.  
  90.     //    Save output at the current point.
  91.  
  92.     if (!exist (output))
  93.     {
  94.       yout.[I] = [t, y.'];
  95.     else
  96.       yout.[I] = output( t, y ).';
  97.     }
  98.     I++;
  99.  
  100.     //    Determine next step size.
  101.     if ( errmax > 6e-4 ) 
  102.     {
  103.       h = 0.9*h*errmax.^-0.2;
  104.     else
  105.       h = 4*h;
  106.     }
  107.   }
  108.   return collapse (yout);
  109. };
  110.  
  111. rk4 = function(ydot, t, h, y, dydt) 
  112. {
  113.   local( k1, k2, k3, k4 );
  114.  
  115.   k1 = h * dydt;
  116.   k2 = h * ydot( t+h/2, y+k1/2 );
  117.   k3 = h * ydot( t+h/2, y+k2/2 );
  118.   k4 = h * ydot( t+h, y+k3 );
  119.  
  120.   return y + ( k1 + 2.0*k2 + 2.0*k3 + k4 ) / 6.0;
  121. };
  122.  
  123. //
  124. // Collapse a LIST of matrices into a single matrix.
  125. //
  126.  
  127. collapse = function ( LIST )
  128. {
  129.   local (i, m, row, col);
  130.  
  131.   row = size (LIST.[members (LIST)[1]])[1];
  132.   col = size (LIST.[members (LIST)[1]])[2];
  133.   m = zeros (length (LIST)*row, col);
  134.  
  135.   for (i in 1:length (LIST))
  136.   {
  137.     m[i;] = LIST.[i];
  138.   }
  139.   return m;
  140. };
  141.